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Introduction 


During  the  execution  of  contract  1962B-7B-C0202*  ARCON 
staff  worked  with  and  provided  computer  software  and  analysis 
support  for  RADC  technical  personnel  on  numerous  projects. 

Most  of  the  analysis  and  software  development  work  involved 
the  following  topics:  circular  antenna  array  theory;  detection 
and  tracking  of  low  flying  targets  by  ground  based  radars; 
null-filter  mobile  radar  software  for  a  CSP-30  system*  detection 
of  a  constant  target  in  log-normal  clutter;  finite  antenna  su¬ 
barrays  irradiated  with  nonorthog ona 1  illumination;  calculations 
of  field  and  radiation  power  patterns  for  a  dipole  array  embedded 
in  a  wire  mesh;  radar  data  processing*  simulation  library  for 
ground-based  radars;  analysis  of  single  wire  individual  resource 
protection  (UIRP)  sensors;  computation  of  received  power  from 
OTH  radars;  LORAN-r elated  programs  and  problems;  multi-hop  pro¬ 
pagation  studies;  modeling  of  VLF/LF  propagations;  whispering 
gallery  theory;  magnetostat ic  surface  wave  analysis*  computa¬ 
tions*  and  plotting*  3-D  plotting  VLF/LF  ionosounding  data; 
ducting  in  the  ionosphere;  ionogram  construction;  magnetic  di¬ 
pole  radiation  caiculat ions*  the  analysis  of  bulk  wave  devices; 
and  an  assortment  of  related  problems.  Representative  examples 
of  the  tasks  performed  will  now  be  presented. 
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■*'j FT. 


1.  Sky-Wav*  Signal  Analysis  and  Calculations 


1.  1  Introduction 

Experimental  methods  for  sounding  the  ionosphere  have  been 
described  in  considerable  detail  elsewhere.  Ref.  1  provides  an 
extensive  list  of  such  references.  Of  interest  to  us  uias  the  de¬ 
termination  of  the  sky-wave  return  signal  based  on  assumed  models 
of  the  ionosphere  and  specification  of  a  transmitted  pulse. 
Because  the  models  are  often  specified  in  the  frequency  domain* 
to  obtain  the  sky-wave  signal  one  is  required  to  compute  Fourier 
Transforms.  Some  of  the  special  issues  that  arise  with  regard  to 
reflection  coefficient  and  full-wave  models  had  to  be  taken  into 
account  as  will  become  evident.  Because  of  the  amount  of  time 
required  to  perform  the  Fourier  transforms  using  the  classical 
computation  approach*  recourse  to  the  Fast  Fourier  transform  al¬ 
gorithm  was  deemed  desirable. 

1.2  Statement  of  the  Problem 

Under  sufficiently  quiescent  conditions*  it  is  reasonable  to 
assume  that  the  ionosphere  does  not  significantly  change  its  fre¬ 
quency-dependent  propagation  characteristics  during  one  or  more 
transmission-reception  time  frames.  Apart  from  additive  noise* 
one  is  justified  in  linearly  relating  the  transmitted  pulse* 
C(t)*  to  the  sky-wave  return  signal*  o(t)*  through  the  convolu¬ 
tion  relationship 

on >  =  at)  @  Utt)  <» 


h(t)  corresponding  to  the  impulse  response  of  the  ionosphere  dur¬ 
ing  the  locally  quiet  period  and Q)  denot ing  aperiodic  convolu- 
t  i  on. 
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With  regard  to  ionospheric  modelling#  the  evaluation  of  (1> 
is  frequently  of  interest.  We  discuss  the  implementation  of  an 
off-the-shelf  Fast  Fourier  transform  (FFT)  code  to  obtain  o(t># 
and  provide  results  based  on  an  analytical  expression  for  the 
frequency  response  of  a  slab  ionosphere  as  well  as  data  generated 
using  a  full-wave  model. 


Equation  (1)  can  be  expressed  in  the  frequency  domain  as 

oO 


Ott) 


=  \  eT^dcJ 


(2) 


-oO 


where  C  and  H  are  Fourier  transforms  of  c(t)  and  h(t>. 

It  is  useful  to  partition  the  integration  interval  (  —  oOjoO  ) 
into  3  regions: 


where 

C  )=  Cuu)  \Aucj)&'°* 


*  A 

If  the  product  C(itO)H(i4J>  vanishes  outside  some  interval 
(  —  Oc  *  )  the  outer  two  integrals  are  lero.  For  the 


t 


\ 

i 


i 


i 

i 


i 
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time-limited  functions  c(t)  and  h<t>*  this  is  strictlg  not  possi¬ 
ble*  but  often  the  approximation  is  verg  good.  The  prolate 
spheroidal  wavef unc t ions  can  be  used  as  a  basis  if  one  is  inter¬ 
ested  in  confining*  in  some  sense  as  "best  as  possible"*  func¬ 
tions  both  in  the  time  and  frequencg  domains  to  finite  bands 

CRef.  23.  But  this  was  not  considered  necessarg  for  the  task 

a  a 

requirements  because  the  product  C<i«)H(iu)  was  nearlg  zero  out¬ 
side  (  -Dc,  Dc  >.  The  central  problem  was  to  develop  a  code 
for  expeditiouslg  computing 

\  A  A  -  ( i 

o  (*)  =  \  c  (CD)  w  (c o)  e  du)  (4) 


A 

for  various  models  of  H(i*>>*  the  ionospheric  transfer  or  frequen¬ 
cg  response  function. 

1.  3  Numerical  Fourier  Transforms 

A 

Because  H<i*>)  is  often  available  as  equispaced  data  samples* 
<4)  was  approximated  with  a  weighted  sum*  i.  e.  *  quadrature  ex¬ 
pression.  It  should  be  noted  that  C <  i*>  )H<  i <•>  )  is  an  aperiodic 
function  in  general.  The  standard  procedure  for  constructing  an 
equispaced  quadrature  approximation  to  (4)  involving  an  aperiodic 
function  is  to  sample  uniformlg: 

=  O*  s  f\  c*=  o>  ±if  ±if ... }  ±  x:  )  (5) 
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It  was  further  assumed  that  o(t)  would  be  evaluated  at  uniformly 
spaced  time  intervals  »  r» *  It  then  follows  that 
with  quadrature  weights  *  * 


f  *  A  -ilSt* 

j  C(itS)  M  (c'o)  €  do  (6) 


* 

C(tu~)  M u'^e. 


n=-£ 


or  in  more  convenient  notation 


O*,  =  Atd  /  %  C„  M„e 


21 fA. 


n=-  *Vi 

The  angular  sampling  frequency  ift  now  given  by  45scxTr/oTt  *  and 
must  be  consistent  with  the  sampling  theorems'  requirement  that 
CJc  £  Ck  f*Z  ■  For  computational  purposes*  one  can  always  choose  a 
higher  cut-off  frequency  lOg*  #  than  the  actual  A)c  associated 
with  the  data*  by  filling  samples  in  the  range  with 

zeros.  This  is  done  so  that 

C3c.  At  ^  t'  <»» 


Under  these  conditions*  and  setting 


6.  =  %  c.  M. 


the  quantity  to  be  computed  takes  the  form 


On  e*fl- 


TTTt  MA 


(H) 


Several  problems  arise  using  (11).  It  is  easy  to  verify 
that  O, r  for  Jt*c*  *  1^1*,  •-  -  By  making  use  of  discrete 
orthogonal ity i  (11)  can  be  inverted  and  also  shows  the  periodic 
property  that  for  The  original 
aperiodic  functions  o(t)  and  £/*!&>>  HfiV*)  have  now  been  trans¬ 
formed  into  periodic  functions  by  virtue  of  the  complex  exponen¬ 
tials  and  equispaced  sampling  in  both  domains.  Furthermore# 
whereas  the  original  aperiodic  functions  correspond  to  different 
frequency  or  time  samples  throughout  their  ranges#  the  induced 

A 

periodicity  of  and  O*  creates  a  discontinuity  at  their 
endpoints  i*  -  z  n/i  the  points  do  not  cor¬ 
respond  to  distinct  times  or  frequencies  because 


rrrc  m 


(±  »/*) 


(12) 


is  independent  of  sign.  Fig  1  depicts  the  situation.  Because 
of  the  periodic  replication  of  the  original  function#  there  will 
be  discontinuities  at 

rv\  («*■  *\\  r  i  ^  =  Oj —  ) 

In  general#  then#  there  can  be  discontinuities  in  the  function 
values  and  derivatives  at  the  endpoints.  With  a  function  discon¬ 
tinuity  in  Ofiu)  #  there  will  be  a  corresponding  component  in 
the  time  domain  function  that  will  fall  off  very  slowly)  namely# 
as  1/m.  To  remedy  this#  one  can  initially  subtract  off  a  linear 
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component  from  Of*u)  end  treat  its  finite  Fourier  transform  sep 
arately  and  analytically  CRef.  33  The  remainder  function  will 
then  fall  off  at  least  as  fast  as  1/m  .  and  can  be  treated  numer¬ 
ically  using  an  FFT  implementation  of  a  DFT.  Other  methods  are 
also  available  to  handle  endpoint  corrections  for  step,  slope  or 
higher  order  discontinuites  CRef.  3  3.  They  were  judged  unneces¬ 
sary  in  our  efforts  because  it  was  found  that  frr  the  data  being 
considered 


j  I  o  L«e?  ° 


To  avoid  counting  the  endpoint  contribution  at  N/2  twice, 
resulting  from  the  periodic  replication,  one  can  set  =°  or 
assign  a  weight  of  1/2  to  both  Ow/x  ant*  ^ee 
Harris  CRef.  43.  For  the  former,  this  leads  to  the  discrete 

Fourier  transform  (DFT)  expression 


cu 


*5  -w/x 


o. 


exp 


£  -  xtt  i  r~\  * 


(13) 


It  is  often  more  convenient  to  make  use  of  the  periodicity 
and  use  a  different  range  of  m  and  n;  viz.  » 

* 

exp  £-1  fft  m<»  /m  J  (14) 

As  shown  in  Fig.  2  for  the  interval  3  .  the  range 
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tM-0  corresponds  to  Cd  “t  #  while  the  range 

corresponds  to  cj  *r  ■£  An  FFT  is  particularly 

convenient  to  use  for  the  determinat ion  of  when  N  >  128. 

Using  1024  complex  samples  in  much  of  our  work#  the  time  required 

to  compute  an  FFT  was  found  to  be  less  than  O.  15  seconds  using 

the  CDC-6600  computer. 

In  the  numerical  Fourier  transform  program  that  was  tailored 
to  sky-wave  data#  most  of  the  time  required  to  obtain  went 

into  "overhead'*  operations:  Construction  of  the  functions  that 
comprise  #  choosing  particular  program  options#  putting  the 

data  into  arrays  of  special  lengths#  windowing  various  functions# 
etc.  Under  these  conditions#  to  obtain  o*,  using  a  string  of 

A 

full-wave  data  for  H(i<J)  and  a  given  c(t)  pulse  took  about  6  sec¬ 
onds  of  CP  time.  By  way  of  comparison#  it  is  interesting  to  note 
that  direct  evaluation  at  only  13  time  points  of  h(t)  for  2000 

A 

samples  of  H(itJ  >  required  about  26  seconds  of  execution  time# 

A 

whereas  using  an  FFT  on  2048  samples  of  H(i*J)  gave  all  2048  sam¬ 
ples  of  h(t)  while  requiring  about  2.3  seconds  of  execution  time. 

1.4  Fast  Convolution 

A  fast  convolution  code  for  continuous-time  functions  was 
written  and  tested  on  cases  that  could  be  evaluated  in  closed 
form.  A  standard  method  for  effecting  the  convolution  was  adopt¬ 
ed  and  will  now  be  summarized.  Variations  of  it  appear  in  the 
literature  and  depend  primarily  on  accuracy  and  speed  require¬ 
ments#  as  well  as  the  class  of  functions  involved.  The  procedure 
we  adopted  was  the  following: 

1.  The  sampled  input  time  signal#  C*  #  is  quadrature 
weighted  using  the  factors  . 

2.  The  DFT  of  ^C,  is  computed  using  an  FFT  code. 

3.  The  frequency  response  or  transfer  function  of  the  sys- 
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tern  (i.e.  #  ionosphere)  is  specified  either  from  a  tape 
file  or  formulas  evaluated  within  the  program. 

4.  The  product  of  the  transformed  input  signal  and  frequen¬ 
cy  response  is  computed. 

5.  The  product  is  smoothed  using  Lanczos  CT-factors*  and 
the  inverse  DFT  is  computed  to  give  the  output  signal. 


The  achieve  aperiodic  convolution*  the  duration  of  the  C<t) 
pulse  and  the  effective  duration  of  the  impulse  response  of  the 
system  must  be  less  than  1/2  the  number  of  DFT  sample  points 
times  the  sampling  interval*  If  this  condition  is  not  met* 

the  output  signal  is  not  convolved  aper iod ical ly  C See  Oppenheim 
and  Schaffer*  Ref.  53.  To  insure  that  this  mas  met*  the  impulse 
response*  h(t)*  was  usually  computed  first  with  a  numerical  Four¬ 
ier  transform  code  and  examined.  This  provided  us  with  informa¬ 
tion  relevant  to  the  effective  time  duration  of  the  impulse  res¬ 
ponse  function. 

For  discrete  functions  the  OT-fac  tors  take  the  form 

/(tV 

O?  -  OJ.,  -o. 

%  -  l. 


*  *»*»*•••,  - 


*-» 


) 


(15) 
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remembering  that  for  an  L-point  DFT#  the  samples  indexed  by  L/2  < 
i  <  L  corresponds  to  the  negative  frequency  (or  time)  region 
The  use  of  0"-factors  is  usually  employed  for  finite  Fourier  ser¬ 
ies.  for  discrete  or  continuous  functions#  especially  if  the 
transformed  function  is  suspected  of  exhibiting  strong  Gibb's 
phenomena#  which  was  the  case  when  a  slab  ionosphere  model  was 
used.  See  Appendix  I.  The  impulse  response  for  the  ionosphere 
in  this  situation  was  observed  to  have  a  “switch-like*1  behavior# 
being  “switched  on”  for  about  50  microseconds  and  off  thereafter. 
The  pulse  or  switch  behavior  was  quite  noticable  for  normal  inci¬ 
dence  to  the  slab  ionosphere.  Using  the  C  -factors#  severe  ring¬ 
ing  effects  were  suppressed#  accompanied  by  negligible  pulse  bro¬ 
adening. 

As  an  alternative  means  for  suppressing  ringing  and  ripples# 
a  4-term  Blackman-Harr is  window  was  also  tried.  Uhile  it  greatly 
suppressed  the  out-of-band  ripples  in  the  pulse-like  h(t)  func¬ 
tion/  it  was  at  the  expense  of  considerable  broadening  and  dis¬ 
torting  of  the  pulse.  The  Blac kman-Harr i s  window  coefficients 
had  been  obtained  from  a  paper  by  Harris  CRef.  43#  and  the 
weights  were  good  to  only  5  significant  digits.  This  created  ac¬ 
curacy  problems  when  the  number  of  DFT  samples  was  increased 
above  1024.  In  the  end#  this  particular  filter  was  eliminated 
from  out  work. 

Because  of  the  nonrecursive  filtering  methods  adopted  in  our 
investigations#  the  broadening  of  the  pulse#  h(t)#  makes  it  ap¬ 
pear  as  if  the  ionosphere  begins  to  respond  prior  to  t  -  0. 

This,  of  course#  is  not  the  case#  rather#  the  non-causal  appear¬ 
ance  is  due  to  the  limitations  of  the  mathematics  in  the  vicinity 
of  the  sudden  response  at  ts0. 
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APPENDIX  I 

Slab  Ionosphere  Model 


For  the  slab  ionosphere#  the  f requenc y-dep endent  reflectivi¬ 
ty  for  TM  radiation  is  of  the  folleming  form: 


Q(iti)  H.  i  — 


l  2d  © 


J 


l  —  ^ 


t  X  KL  Gn& 


(A-l) 


where 


Q(<v)  - 


K.  - 


-L1+  issi  1 4  fl+  jfcj 

icA'+  *#cr/ c^ee  1 


(A- 2) 


XTT 


-f/c  r  /0/C 


^3"  r  C#«  Jwt+ivi+y 


*•5-  E-^ 


ih!\aS  J** 


L  ■  $Ut  ‘H*t'cfc*#ss  e-^.  7-^ 

€©  =  I-  €  -  *1  /at-* 
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To  avoid  the  "apparent"  singularity  at  4)-  0#  R  (0)  was  evaluated 
numerically  at  4)  *  0*  s  1.  E-9.  Using  Rm  as  the 
frequency-response  of  the  ionosphere*  the  impulse  response  was 
then  determined  through  a  numerical  Fourier  transform  computa¬ 
tion. 
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2.  Bulk  Wave  Investigations 


2.  1  SSBW  Transducers 

Surface-skimming  bulk  wave  (SSBW)  transducers  were  of  inter¬ 
est  to  us  in  the  context  of  analysis  and  design.  An  interest  in 
materials  and  crystal  cuts  that  would  support  certain  acoustic 
modes  for  filter  applications  led  us  to  solve  the  Christoffel 
equations  for  the  allowed  propagation  c  ?**»“*-  *ts  similar  to  the 
analyses  developed  by  Milsom  CRef.  \ Hussein  and  Ristic 
CRef.  23.  The  condition  for  allowed  ^  ;s  that  the  determi¬ 
nant  of  the  Christoffel  equations  vam a  given  surface  wave 
velocity*  With  this  information  could  proceed  to  satis¬ 
fy  the  inhomogeneous  boundary  conditions*  in  k-space<  for  a  line 
source.  Our  intermediate  objective  w*i’^  to  be  able  to  arrive  at 
the  effective  permittivity  function  for  an  material  and  crystal 
orientation  specified. 

Because  of  the  goal  to  treat  arbitrary  crystal  orientations* 
it  was  necessary  to  develop  or  obtain  software  that  would  be  gen¬ 
eral  enough  to  express  the  4th*  3rd*  and  2nd  rank  tensors  (stiff¬ 
ness*  piezoelectric*  and  dielectric  tensors*  respectively)*  for 
any  particular  crystal  cut*  in  the  Christoffel  and  boundary  con¬ 
ditions  equations.  An  existing  program  called  LAYERS  CRef.  33 
appeared  to  fill  many  of  these  needs*  and  it  was  decided  to  make 
use  of  as  much  of  it  as  we  could. 

Our  efforts  were  stymied  very  early  using  the  LAYERS  code 
because  the  Christoffel  equations  did  not  yield  slowness  surfaces 
in  agreement  with  published  curves.  Qual i tat i vel y  they  were  sim¬ 
ilar*  but  the  magnitude  of  the  deviations  were  intolerable  when 
compared  with  the  work  of  Milsom*  Wagers*  and  Ristic  CRef.  2* 
43.  There  was  no  way  to  scale  our  results  to  agree  with  theirs. 

Some  of  our  slowness  curves  appeared  to  exhibit  symptoms  of 
numerical  instability  and  poor  convergence.  Tolerances  were 
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tightened  and  loosened  to  very  little  avail.  We  then  examined 
the  root-finding  code  (Muller's  inverse  parabolic  Interpolation 
algorithm)  for  evaluating  the  determinant#  and  convinced  our¬ 
selves  that  it  was  among  the  very  best  off-the-shelf  codes  one 
could  use.  It  converged  rapidly  and  accurately  in  actual  tests. 
The  Christoffel  matrix  itself  was  examined  as  best  we  could#  con¬ 
sidering  that  rotated  tensor  expressions  were  invjlved.  Some  un¬ 
usual  scaling  was  noted  in  the  code  but  we  wer^  unable  to  ascei — 
tain  its  possible  role  in  generating  errors.  The  project  was 
eventually  terminated. 

2. 2  Power  Patterns  for  SH-Wave  IDTs 

2.  2.  1  Analysis 

Because  of  the  potential  applications  of  microwave  frequency 
array— steer ed  bulk  wave  devices#  it  was  decided  that  the  radia¬ 
tion  characteristics  of  shear  horizontal  (SH)  waves  should  be  in¬ 
vestigated.  One  quantity  of  particular  interest  that  needed  to 
be  computed  was  the  far-field  power  pattern  of  a  planar  array  of 
electrodes.  The  derivation  of  it  is  carried  out  below.  Fig. 
3  shows#  for  a  double-electrode  array  and  relevant  geometry. 
As  several  investigators  have  shown#  radiation  into  an  isotropic 
medium  is  often  a  reasonable  physical  approximation  to  make*  and 
certainly  makes  the  analysis  simpler. 

For  the  double-electrode  array  shown  in  Fig.  3  *  we  note 

that 


Const. 
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Expanding  {  leads  to  the  following: 

^  0  -  Z.Rx,>4A0  +  xM 

r  z /u^e  4  -  ■£*  h*«*0 

+  (x*w  +  *T  ~  *£fc  ****♦»)  i 


For  the  far-field  power  distribution#  i.e. 


U  *  - 


«  i  > 


the  quadratic  term  is  negligible  by  comparison  It  follows  that 


f  =T^  j  I  -a**0  (^w  +  i  -k  “  ®*WO+*-  \ 


By  far-field  we  also  mean  that  the  half-aperture  of  the  transmit¬ 
ting  array  is  much  less  than  the  distance  to  the  observation 
point  at  R. 


1+  ■\r?»*T3 


Following  Farnell's  analysis  CRef  1*  it  is  clear  that 
the  far-field  potential  due  to  all  the  line  sources  comprising 
the  array  is  given  by 


T/X 


(6) 


i  r  -3/x. 


where 
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Using  the  asymptotic  form  for  the  Hankel  function  gives 
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Denoting  the  integral  by 

=  jd<-  'T.,»)glk,iw~'e 

19 


I 


and  expanding  the  charge  density  as  usual  CRef.  1 »  i.  e.  * 

~/L  CT  X<^)  <“» 

t\~  O 


leads  to 


=  -rr  "21  Cm**  (-i)  X  (?)  (H) 

^  v\*0 


where 


<3”*  =  4=/ia 

and 

=  -sA  • 

2. 22  Results  and  Discussion 

For  J  =  200#  i.  e.  #  100  pairs  of  electrodes#  Fig.  shows 

the  far-field  power  pattern  for  =  1/2.  For  an  isotropic  sub- 
strate  y'fV)  =  ¥  (-1)  This  is  not  the  case  for  an 

anisotropic  medium#  but  the  difference  has  been  shown  to  be  quite 
small  in  many  situations  by  Lewis  [Ref.  51.  It  is  clear  that 
the  DE  array  tuned  to  its  fundamental  frequency  behaves  like  an 
end-five  antenna.  For  these  computations  the  array  was  assumed 
to  have  a  periodic  polarity  impressed  on  the  pairs  of  electrodes. 
The  infinite  array  approximation  was  used  in  the  sense  that  the 
expansion  coefficients  were  the  same  (apart  from  a  sign)  for 
all  electrode  pairs.  The  expansions  were  calculated  using  N  =  6# 
and  the  coefficients  had  the  following  magnitudes: 

C(t  *  0.  1453 14 
=  0.  059185 
=  0.  002699 
Cj  =  o.  001402 
=  0.  000047 
=  0.  000031 
=  O.  000001 

Element  factors  based  on  calculations  and  derivations  developed 
above  will  be  determined  for  both  simple  and  double-electrodes# 
and  will  appear  in  a  journal  publication. 
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Far- Field  Power  Pattern 
100  Electrode  Pairs 

Figure  4 
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3.  Message  I  dent i f i cat i on 


A  low-f requenc y  method  of  communication  was  tested  between 
Forestport*  New  York  and  Troy*  New  Hampshire.  This  method  in¬ 
volved  the  transmission*  reception  and  subsequent  analysis  of 
time-delayed  low-frequency  groundwaves  <25  KHz)  which  were  digi¬ 
tized  and  stored  on  magnetic  tape.  The  analysis  detailed  below* 
included  determining  the  time-delay  of  the  groundwave  in  each  re¬ 
cord.  A  message  was  encoded  by  varying  the  time-delay  of  a 
number  of  consecutive  records  wherein  each  delay  corresponded  to 
certain  alphabetic  character. 

A  computer  program  was  designed  to  perform  correlations 
between  a  standard  waveform  and  that  of  each  record  on  tape  by 
mathematically  shifting  the  former  along  the  latter  to  determine 
the  position  of  maximum  correlation.  This  point  represents  the 
time-delay  of  the  given  waveform. 

A  well-known  cross-correlation  function 


X  x*  *+  ?  Z.  *<-  +  *  Z. 

X(K\  = 

[_  2T  K  -  Z  -  *  ** 


+  *  Z  4  -  (*♦*> 

q’A. 

%  -  *  V*  J 


was  used  in  which  x  and  y  are  the  values  of  the  standard  waveform 
and  the  one  being  tested*  n*  is  the  total  number  of  elements  in 
each  waveform  and  k  is  the  particular  shift  of  the  standard  along 
the  waveform  being  compared. 

The  field  tapes  consisted  of  records  of  length  256*  each 
representing  a  time-frame  of  the  same  number  of  microseconds  and 
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each  containing  one  waveform  40  microseconds  long.  Blocks  of  22 
records  were  created  bg  transmitting  22  consecutive  waveforms 
with  the  same  time-delag.  Mang  tapes  containing  waveforms  of 
varging  signal  strengths  were  created  for  later  analgsis. 

Initial  results*  comparing  the  standard  waveform  to  those  of 
a  tape  of  high  signal-to-noise  ratio  were  successful*  i.  e.  *  the 
message  was  easilg  decoded.  There  were  problems*  however*  decod¬ 
ing  tapes  having  low  S/N  ratios.  Averaging  five  consecutive  re¬ 
cords  before  cross-correlating  was  tried  but  since  there  is  no 
physical  indicator  of  the  beginning  of  a  given  block*  it  was  im¬ 
possible  to  be  sure  whether  one  was  averaging  records  in  the  same 
block  (same  time-delags)  or  those  in  two  different  blocks. 
Unders tandab ly,  averaging  records  in  two  different  blocks  would 
tend  to  worsen  results  since  an  additional  waveform  would  be  in¬ 
troduced  into  the  averaged  record  and  this  was  found  to  be  the 
case. 


It  was  determined  that  more  advanced  signal  processing  tech¬ 
niques  are  necessarg  to  improve  the  S/N  ratio  prior  to  applging 
the  cross-correlation  method.  Current  investigations  include  the 
use  of  fast  Fourier  transforms  and  filtering  to  eliminate  noise 
and  improve  the  S/N  ratio. 
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4.  Nonlinear  Transcendents  1  Equations 


Our  objective  was  to  solve  a  system  of  5  nonlinear  transcen¬ 
dental  equations.  The  first  3  pertained  to  the  reflection  of  an 
incident  ray  from  a  flat  facet  having  a  slope  angle  ^  and  an  ori¬ 
entation  angle  at  with  respect  to  the  transmi tter-rec ei ver  line  of 
sight.  The  last  two  equations  ensured  that  the  reflected  ray  ar¬ 
rived  at  the  receiver  if  the  height  of  the  ^oint  of  reflection 
was  required  to  be  H  and  the  distance  between  the  transmitter  and 
the  receiver*  D. 

An  attempt  to  solve  the  equations  analytically  failed.  The 
next  effort  was  to  obtain  an  approximate  solution  of  the  system 
of  equations.  The  system  was  transformed  through  several  substi¬ 
tutions  into  a  system  of  algebraic  equations  with  the  hope  of  di¬ 
minishing  the  computer  time  demands#  because  it  was  known  that 
searching  for  the  solution  using  an  iterative  process  would  be 
c omputat ional 1 y  expensive#  especially  when  there  was  no  informa¬ 
tion  for  good  starting  values.  A  rapidly  converging  iterative 
process  was  applied  to  the  algebraic  system#  the  subroutine  ZCNT 
from  the  IMSL  library  was  used.  However#  it  turned  out  to  take 
considerable  computer  time#  2000  sec.  For  this  reason#  the 
paper  from  which  these  equations  were  derived  was  reexamined  and 
found  to  have  a  misprint.  A  careful  examination  of  the  geometry 
showed  the  equations  weren't  necessary  at  all  for  the  determina¬ 
tion  of  the  directions  i  ^,6  >  of  the  incident  ray  and  direc¬ 
tions  (  >  of  the  reflected  ray.  All  unknown  angles  could  be 
found  directly  from  the  geometry  of  the  reflection  by  solving 
simple  solid  geomtery  problems. 

The  new  calculations  took  less  than  1  second  of  computer 
time#  and  the  results  had  obvious  geometrical  meaning.  Plots 
were  made  of  for  different  values  of  the  parameters. 
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5.  Tropospheric  Range  Error  Correction  Utilizing 
Ex oatmospher i c  Sources 

5. 1  Introduction 

Tropospheric  refractive  range  and  angle  errors  limit  the 
performance  of  radar  and  navigation  systems  that  operate  at  low 
elevation  angles.  For  some  locations  and  for  many  applications# 
antenna  pointing  corrections  based  on  surface  refractivity  alone 
or,  if  necessary#  a  vertical  refractivity  profile  are  adequate. 
However#  new  systems  that  operate  at  very  low  elevation  angles 
r^q^uire  improved  accuracy.  More  accurate  angle  error  corrections 
can  be  obtained  by  utilizing  calibration  sources  (such  as  the 
limb  of  the  sun#  radio  sources#  satellites)  the  angular  positions 
of  which  are  precisely  known  to  a  sufficient  accuracy — the  range 
error  correction  is  still  a  problem.  It  has  been  shown  that  it 
is  possible  to  develop  a  method  for  obtaining  an  improved  range 
error  correction  by  utilizing  measured  angle  error  data  for  the 
same  ray  path.  An  expression  for  the  refractive  angle  correction 
for  a  target#  based  on  the  measured  angle  error  of  the  calibra¬ 
tion  source  is  first  derived  <1>. 

Both  the  target  and  calibration  source  are  assumed  to  be 
outside  the  troposphere  where  the  index  of  refraction  n(r)  is  ap¬ 
proximately  unity#  hence  the  angular  refraction  T  is  the  same  for 
both 

/ u w  ( r -  £>)  L  m+MW  c«*V0a-€,>  -  a.  a^(0+-€. 

(1) 

r  ^  (T-£x\\ 
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In  this  equation,  #  the  angle  error  of  the  calibration  source#  ! 

and  the  initial  elevation  angle  of  antenna  are  known*  as  are  ap¬ 
proximate  heights  h,  and  h^  of  the  target  and  calibration  source 
(see  Fig.  5  ). 

i 

There  is  a  direct  formula  to  express  the  angle  error  correc-  I 

tions  €|  #  (for  the  target  or  for  the  calibration  source  res¬ 
pectively)  through  the  angular  refraction  j 


j  Ga/T-  aL*'  1 

[  <u$>  /n*  •  ”tL  —£^&t  J 

(2) 

where 


An  J 


-  » 

Bi  - 


r. 

i 


under  assumptions  mentioned  above. 
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However#  before  th 0  angle  error  of  the  target  can  be  de¬ 
termined*  the  angular  refraction  X  must  be  found 


-  A i~-  L«.*kT 


(3) 


The  angular  position  of  calibration  sources  are  precisely 
known*  so  it  is  possible  to  calculate  the  angle  correction  Cx  by 
taking  the  difference  of  the  apparent  position  which  is  meas¬ 
ured  and  the  true  position  which  is  tabulated: 

£*  -  Q«  .  (4) 


♦It  should  be  emphasized  that  the  following  formula  is  correct  as 
well. 

X  =  -"Al  4-  -  CnT’b^ 

(31) 

-  ^-’LcSC  Ca^-^.Q 


since  the  angular  refraction  X  for  the  target  is  the  same.  Then 
with  X  known  from  <3>  the  angle  error  of  the  target  €,  can  be 
calculated  from  <1>. 

30 

-  ■  - - - 

**  I*  *  Jr.r  *  •  v 


It  is  then  shown  that  the  range  error  correction  for  a  ar- 
get  can  be  expressed  as  a  function  (integral)  of  the  angular  re¬ 
fraction.  The  angle  error  for  the  target  is  also  involved. 
Therefore*  the  range  error  correction  can  be  calculated  from  a 
set  of  angle  error  measurements  (4). 

The  range  error  AR*  <a#  Cfc#  rg  )  correction  is  defined  bg 


rt)  =  ~Re(4,C«,r,)  (5) 


where  Rcand  Re  are  the  apparent  (or  radio)  range  and  the  true  (or 
slant)  range  respectively. 


The  expression  for  calculating  the  range  error  correction  is 
given  by 


C 


with  following  denotations: 

-  is  where  Q  -  the  elevation  angle  of  the 

antenna; 

C  -  is  cos  <0)»  where  $  the  elevation  angle  of  the 
calibration  source; 

-  the  distance  of  the  antenna  from  the  center  of 

the  earth; 

*  T,  -  the  distance  of  the  target  from  the  center  of  the 
earth; 

fi*  -  the  quotient  (a  +  h(  )/a; 

-  the  angle  error  of  the  target; 

Hi  (qjC+j  r,  )  “  is  the  total  angular  refraction  along  the  ray 
path  (the  ray  path  lies  in  the  great  circle 
plane#  determined  by  the  locations  of  the  anten¬ 
na  and  target#  see  Fig.  1). 


T  r, )  =  - 


The  function  n<r)  is  expressed  as  follows: 


0(0  =  I  +•  KJ(r-A.)  *  jo* , 

r  u,  4  aw  *  (*-«> 
N<*--0  -  \  KJ, 


loS  L~  0.142-4  (r- a- 


o  4  *-«.  i  I 
\  4  r-o>  4 

•V  <  r~*' 


\ 

\ 


32 


1 

I 

!  1 

J 

here  j  ] 

N>$  s  "ill  .  o  i 
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KJ,  =  Kir  ♦  AU 

C  -  ijiu/s) 


This  formula  for  N<r-a)  corresponds  to  "Modified  Effective 
Earth's  Radius  Model"  which  is  usually  used  to  solve  refraction 
prob  1  ems. 

It  should  be  emphasized  that  (6)  provides  the  difference  in 
range  error  for  elevation  angles  corresponding  to  the  limits 
and  C.  For  example*  for  the  case  of  the  sun  as  a  calibration 


source* 


it  is  the  difference  of  range  errors  for 


and 


corresponding  to  elevation  angles  of  60  and  *  respectively. 

In  order  to  obtain  the  total  range  error  for  an  elevation 

angle  C*.  <  $  f 1  )  it  is  necessary  to  estimate  the 

o 

total  range  error  for  an  angle  of  60  .  This  is  not  considered  a 
serious  limitation  since  the  angular  refraction  for  very  high  an¬ 
gles  is  minimal*  also*  the  corresponding  range  error  is  rela¬ 
tively  small. 

5.2  Task  Description 

The  problem  to  analyze  consisted  of  the  following  steps: 

1)  to  select  the  numerical  integration  technique  which  provides 
the  least  number  of  measurements  required  to  guarant  *  a 

given  accuracy  of  integration  for  the  integral 
C* 


\  Z  ( a,C*,  r, )  <IC* 


(8) 
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2)  to  analyze  the  sensitivity  of  this  integral  to  the  small  pei — 
turbation  of  Q  values  (the  measurement  errors). 

3)  to  analyze  the  sensitivity  of  the  angle  error  €,  to  the  small 
perturbations  of  0  values  (the  measurement  errors). 


4)  to  analyze  the  sensitivity  of  the  range  error  correction 
A |  to  the  small  perturbat ions  of  ^  values.  Values 

of  in  expression  (3)  are  supposed  to  be  exact#  so 


A€ ,  =  aQc 


To  obtain  the  range  error  correction  from  (6)  the  following 

procedure  was  suggesed: 

1)  Select  a  set  of  elevation  angles  between  and  0  and  obtain 
the  correspond ing  set  of  the  angle  errors  of  the  calibration 
source#  utilizing  the  relation  (4).  Elevation  angles  Q  are 
to  be  measured#  the  appropriate  true  angular  positions  p  of 
the  calibration  source  are  to  be  found  from  tables  <if^  for 
example#  the  sun  is  the  calibration  source#  the  solar  ephem- 
eris  is  to  be  utilized). 

2)  Calculate  the  appropriate  angular  refractions  for  the  se¬ 
lected  set  of  elevation  angles  from  the  relation  (3). 

3)  Calculate  the  angle  error  £,  from  the  equation  (1).  For  this 
calculation  only  one  measurement  is  needed. 

4)  Use  the  obtained  set  of  X  values  for  the  approximate  numeri¬ 
cal  calculation  of  the  integral  on  the  right  side  of  the  ex¬ 
pression  (6). 
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The  accuracy  of  this  technique  is  dependent  on  many  factors 
Especially#  of  primary  importance  are  the  accuracy  with  which  the 
angle  error  can  be  measured#  the  refractivity  structure  of  the 
troposphere  and  the  number  of  angle  error  measurements. 

Of  secondary  importance  are  the  accuracy  of  the  heights  and 
angular  coordinates  of  the  calibration  sources#  the  height  of  the 
target#  the  measurement  of  the  surface  index  of  refraction  # 

the  numerical  integration  techniques  that  are  employed#  and  the 
estimate  of  the  total  range  error  contribution  from  the  angular 
refraction  at  elevation  angles  above  the  maximum  elevation  angle 
of  calibration  sources. 

5. 3  Work  Description 

1)  A  numerical  experiment  was  designed  to  select  the  appropriate 
numerical  integration  technique  which  would  guarantee  a  given  ac¬ 
curacy  with  the  least  number  of  nodes  #  i.  e.  #  that  would  use 
the  least  number  of  required  measurements.  The  exact  theoretical 
expression  for  the  angular  refractivity 

\  r/t»>  _ 

an*C-J  (91 

a 

was  utilized  in  this  experiment. 

It  should  be  emphasized  that  this  integral  is  close  to  an 
integral  with  a  singularity  when  O 

The  relative  accuracy  of  calculations  considered  to  be  equal 
to  0*5  *  i  there  was  no  reason  to  consider  higher  accura¬ 

cies  since  the  value  a#  the  antenna  distance  from  the  center  of 
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the  earth*  had  no  more  than  3  correct  significant  digits.  Three 
methods  for  numer ical  integration  were  c ompared :  Romberg  rule* 
Simpson  formula  and  Hermite  formula,  the  nodes  being  selected  au¬ 
tomatically.  It  turned  out  that  Simpson's  method  gave  noticable 
better  results;  for  example: 


e 

a 

Number  of  Nodes 

Simpson  Pule 

Romberg  Rule 

3* 

25 

41 

• 

o 

fNj 

II 

€L> 

0 

o 

o 

53 

73 

The  nodes  selected  by  Simpson's  method  were  utilized  for 
calculation  of  the  integral  <8>*  but  the  angular  refra ctivity 
was  computed  from  the  expression  (3)  with  the  appropriate  in¬ 
terpolation  along  the  table  Results  of  two  different  ways 

of  the  calculation  of  the  integral  (8)  coincided  with  3  signifi¬ 
cant  digits  Therefore*  the  set  of  angles  obtained  from  the 
first  preliminary  calculation  utilizing  the  expression  (9)  for 
the  angular  refractivity  can  be  used  as  the  set  of  measurement 
ang  les. 

2)  The  sensitivity  of  the  integral  (8)  to  perturbations  was 

checked  numerically.  Three  significant  digits  of  the  integral 
remained  unchanged  while  perturbations  didn't  exceed 
0*5  IQ  in  absolute  values.  Such  precise  measurement 
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can  be  accomplished  practically#  so  the  sensitivity  is  suitably 

1  OU). 


* 


3)  The  sensitivity  of  angle  error  €,  with  respect  to  perturbation 
of  angle  error  of  the  calibration  source#  (under  assumption 
that  are  exact)#  has  been  checked  numerically.  It 
turned  out  that 


I  I  *  z  I  ScJ 

s. 


when  angles  0^  belong  to  the  interval  0#  f  $  S 


t 

i 
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6  Maximum  Signal  Power  at  a  Critical  Down-Range  Distance 

6. 1  Introduction 

Over-the-hor i ion  (OTH)  radar  signals  are  reflected  from  sev¬ 
eral  layers  of  the  ionosphere  before  illuminating  some  point  of 
the  earth.  The  optimum  frequency  -f  and  launch  angle  p  to  be 

utilized  may  be  predicted  via  ray-tracing  techniques.  The  goal 
of  the  program  PWRMAX  is  to  determine  -f  and  (S  for  maximum  signal 
power  at  a  critical  down-range  rc  Both  maximums  were  needed; 

the  local  maximum  of  power  and  absolute  maximum  of  power. 

In  order  to  reach  this  goal  the  program  PURMAX  combines  a 
high-speed  ray-tracing  computer  code  with  a  power  versus  range  r 
curve-fitting  algorithm. 

6  2  General  Discussion 

The  search  for  the  optimum  frequency  is  performed  over  the 
initial  intervals  of  frequencies  CFj  *  3  and  angles  CB,  * 

Bx  3#  the  intervals  being  determined  from  physical  considera- 
t  i  ons. 

The  ray-tracing  computer  code  allows  one  to  obtain  the  cor¬ 
responding  values  of  power  p  and  range  r  for  each  pair 

($,(!)€  L  F,  ,  F,.  ]  <2>  LB,  3  . 

For  each  fixed  frequency  ^  €  ;W)  and 

variable  j  r  . . .  tM  }  a  table  of  "power  p_  11 

versus  "range  r.  "  may  be  obtained: 


Such  tables  represent  functional  dependencies  of  “power"  versus 
"range"  and  must  be  examined  at  the  beginning. 

The  functional  dependencies  mag  be  described  approx imatelg 
bg  the  quadratic  obtained  using  standard  least  squares.  Denoting 

a 

the  appropriate  polgnomial  bg  p  (r)  *  the  following  approximate 

» 

equalitg  is  correct: 


T?  <  n  s.  T?  (O 


(i) 


The  value  of  power  at  the  point  ^  mag  be  approximately  calculat¬ 
ed  as  ?.  (rt)  : 


T?  to 


(2) 


This  leads  to  a  table  of  the  form 


f, 

fi 

•  •  • 

f. 

— X — 

•  •  • 

i 

P- 

which  mag  be  treated  as  an  approximate  representation  of 
functional  dependence  Pc  (f)  “the  value  of  power  at 


the 

the 
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down -rang*  rc  M  versus  "frequency".  Using  the  functional  form  of 
the  above  table*  the  absolute  maximum  may  be  obtained  with  stan¬ 
dard  methods  (involving  either  simple  selection  of  the  least 
square  curve  fitting  with  further  selection  of  absolute  maximum 
for  a  parabola. 

With  regard  to  the  local  maximum  of  power*  the  least  squares 
curve-fitting  usually  is  not  a  suitable  method  to  investigate 
local  properties  of  fitted  curves.  But  if  certain  conditions 
apply*  the  fitting  curve  may  be  used  for  approximate  calculations 
of  the  local  maximum. 

Functional  dependencies  f/r  #,  l, ...  9*t)  were  assumed  to 

be  quite  smooth  and  to  have  the  single  local  maximum  in  the 
neighborhood  of  the  point  rc  This  assumption  was  justified 

physically  and  with  numerical  experiments. 

If  a  curve  is  shaped  as  assumed*  the  coordinates 

(f-  ,  )  of  the  fitting  parab  ola  local  maximum 

may  be  regarded  as  an  approx imat ion  to  the  coordinates 

( T***  >  )  of  the  local  maximum  of  (*>  .  The  func¬ 

tional  dependency  "frequency"  versus  "range-coordinate  of  local 
maximum"*  may  then  be  approximately  restored  using 

standard  least  square  so  that 


A 

Ft  <*•)  c  F,  (n 

Amx 


where  B|  is  a  polynomial  of  second  degree. 

Finally,  the  frequency  that  provides  the  local  maximum  of 
power  at  the  point  I*  may  be  calculated  approx imately  as 

A 

ferf.  *  fit— .tr‘>  (4) 
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7.  Multi-Hop  Propagation 
7.  1  Introduction 


The  computation  of  the  number  of  hops  versus  launch  angle 
for  geometric  multi-hop  propagation  between  concentric  boundaries 
presents  certain  analytical  problems. 

This  problem  is  mainly  due  to  the  implicit  form  in  which  re¬ 
lations  between  number  of  hops*  launch  angle*  relative  transmis¬ 
sion  delay*  distance  and  height  of  transmission  boundary  are  pre¬ 
sented. 

To  obtain  iso-height  and  transmission-delay  contours#  expli¬ 
cit  formulas  for  number  of  hops  and  launch  angle  of  the  ray  were 
der i ved. 


This  representa t i on  also  obviates  the  separate  treatment  of 
multi-hop  round-th e— wor 1 d  propagation  as  an  additional  case  to 
the  propagation  over  a  certain  distance  along  the  surface  of 
ear th . 

7  2  Analysis  of  Parametric  Equations 

For  mirror  reflections  between  concentric  spherical  ionos¬ 
phere  and  earth  over  a  distance  <d)  along  the  surface  of  earth 
the  following  parametric  equations  apply  for  the  height  of  re¬ 
flection  boundary  <h)  and  relative  transmission-delay  <8  ) 


where  R  is  radius  of  earth#  6  is  launch  angle  of  the  rag# 
number  of  hops  and  c  is  velocity  of  light. 


n 


In  addition#  the  round-the-world  case  is  presented  by  equa¬ 
tions  v 


An  analysis  of  given  parametric  equations  yields  a  unique  set  of 
formulas  for  launch  angle  (0)  and  number  of  hops  (h). 

r\  -  -hr*.  "  3 

(5) 


where 


i 

I 


i 


C.^  ■*  J 

O.R 


(6) 


*5  ' 


(?) 


i 


Thus#  the  round-the-wor 1 d  case  could  be  derived  easily  by  substi¬ 
tuting  distance  d  with  2vHt  and  iso-height  and  transmission  delay 
contours  could  be  easily  obtained. 
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